################################################################################
##
## Purpose: This script creates all figures and tables from SI section 10
##
## Author: James Bisbee (james.h.bisbee@vanderbilt.edu)
##
## Input Files:
##  - ./data/prepped/finalData.RData: Prepped data from 9_DATA_final_build.R
##
## Output Files:
##  - ./output/figures/SI_figure_24.pdf
##  - ./output/figures/SI_figure_25.pdf
##  - ./output/figures/SI_figure_26.pdf
##  - ./output/figures/SI_figure_27.pdf
##  - ./output/figures/SI_figure_28.pdf
##  - ./output/tables/SI_table_11.tex
##
##
## See associated log file for compute environment, package versions, 
##  and date of most recent run.
##
################################################################################
rm(list = ls())
gc()
require(fixest)
require(tidyverse)
require(ggridges)
require(marginaleffects)
require(patchwork)
require(cowplot)
require(quantreg)
require(lfe)

set.seed(123)

# Compute details
print(paste0('Compute environment from ',Sys.Date(),' run by Bisbee'))
if(Sys.info()['sysname'] == 'Windows') {
  ram_size = system("wmic MemoryChip get Capacity", intern = TRUE)[-1]
  model_name = system("wmic cpu get name", intern = TRUE)[2] # nocov
  vendor_id = system("wmic cpu get manufacturer", intern = TRUE)[2] # nocov
  
  print(list(ram = stringr::str_squish(ram_size)[1],
             vendor_id = stringr::str_squish(vendor_id),
             model_name = stringr::str_squish(model_name),
             no_of_cores = parallel::detectCores()))
} else if(Sys.info()['sysname'] == 'Linuxs') {
  splitted <- strsplit(system("ps -C rsession -o %cpu,%mem,pid,cmd", intern = TRUE), " ")
  df <- do.call(rbind, lapply(splitted[-1], 
                              function(x) data.frame(
                                cpu = as.numeric(x[2]),
                                mem = as.numeric(x[4]),
                                pid = as.numeric(x[5]),
                                cmd = paste(x[-c(1:5)], collapse = " "))))
  df
} else {
  cat("If not on Linux or Windows, you'll have to figure out your own solution to seeing the compute environment.")
}

sessionInfo()


load('./data/prepped/finalData.RData')



# SI Figure 24
test <- utterance_level %>%
  select(respondingTo,interruptor,matches('SENT_')) %>%
  select(-matches('SENT_.*lag'),-SENT_error) %>%
  gather(sentiment,value,-respondingTo,-interruptor)

pdf('./output/figures/SI_figure_24.pdf',width = 7,height = 5)
test %>%
  group_by(respondingTo,sentiment,interruptor) %>%
  mutate(id = row_number(),
         nUtt = n()) %>%
  filter(nUtt > 2) %>%
  select(-nUtt) %>%
  spread(interruptor,value,sep = '_') %>%
  drop_na() %>%
  mutate(diff = interruptor_1 - interruptor_0) %>%
  group_by(sentiment) %>%
  summarise(diff = mean(diff,na.rm=T),
            pval = t.test(interruptor_1,interruptor_0)$p.value,
            tstat = t.test(interruptor_1,interruptor_0)$statistic) %>%
  ungroup() %>%
  mutate(se = diff/tstat) %>%
  filter(!grepl('autho|incoh|unsubs|reject',tolower(sentiment))) %>%
  ggplot(aes(x = diff,y = reorder(sentiment,diff))) + 
  geom_point() + 
  geom_errorbarh(aes(xmin = diff - 2*se,xmax = diff + 2*se),height = 0) + 
  geom_vline(xintercept = 0,linetype = 'dashed') + 
  labs(x = 'Difference',
       y = 'Tone',
       title = 'Differences between interruptions and other utterances',
       subtitle = 'T-test difference of means by tone') + 
  theme_bw()
dev.off()


# SI Figure 25
m <- feols(SENT_combAttack ~ respondingTo*factor(interruptor) + log(nchars) + log(nchars_lag) + factor(interrupted) | opensecretsID,
           utterance_level %>%
             filter(all > 30,ind > mind) %>%
             filter(!grepl('FED',opensecretsID)) %>%
             filter(grepl('FED',respondingTo)) %>%
             mutate(respondingTo = ifelse(grepl('YELLEN',respondingTo),'Yellen',
                                          ifelse(grepl('FED',respondingTo),'Male Fed Chairs',NA))) %>%
             mutate(respondingTo = relevel(factor(respondingTo),ref = 'Male Fed Chairs')),
           cluster = 'respondingTo')


toplotCAP <- plot_cap(m,condition = list('respondingTo','interruptor'),draw = F)
toplotCME <- plot_cme(m,variables = 'respondingTo',condition = list('interruptor'),draw = F)

pCAP <- toplotCAP %>%
  as_tibble() %>%
  ggplot(aes(x = respondingTo,y = estimate,color = factor(interruptor),
             shape = factor(interruptor))) + 
  geom_point(position = position_dodge(width = .1),size = 2.5) + 
  geom_errorbar(aes(ymin = conf.low,ymax = conf.high),width = .05,
                position = position_dodge(width = .1)) + 
  scale_color_manual(name = 'Utterance',values = c('grey40','black'),labels = c('Non-Interruption','Interruption')) + 
  scale_shape_discrete(name = 'Utterance',labels = c('Non-Interruption','Interruption')) + 
  labs(x = 'Responding To...',
       y = 'Predicted Aggression',
       title = 'Predicted Values',
       subtitle = 'Aggression by utterance') + 
  theme_bw()


(pCME <- toplotCME %>%
    ggplot(aes(x = factor(interruptor),y = estimate,color = factor(interruptor),shape = factor(interruptor))) + 
    geom_point(size = 2.5) + 
    geom_errorbar(aes(ymin = conf.low,ymax = conf.high),width = .05) + 
    geom_hline(yintercept = 0,linetype = 'dashed') + 
    theme_bw() + 
    scale_color_manual(name = 'Utterance',values = c('grey40','black'),labels = c('Non-Interruption','Interruption')) + 
    scale_shape_discrete(name = 'Utterance',labels = c('Non-Interruption','Interruption')) + 
    scale_x_discrete(name = 'Utterance',labels = c('Non-interruption','Interruption')) + 
    labs(y = 'Coefficient Estimate',
         title = 'Marginal Effects',
         subtitle = 'Yellen vs Fed Chairs'))


pdf('./output/figures/SI_figure_25.pdf',width = 7,height = 4)
pCAP + pCME + plot_layout(guides = 'collect') & theme(legend.position = 'bottom')
dev.off()

# SI Figure 26
# First, fine the quantile that corresponds to the standard OLS
which.min(abs(mean(utterance_level %>%
                     filter(all > 30,ind > mind) %>%
                     filter(!grepl('FED',opensecretsID)) %>%
                     filter(grepl('FED',respondingTo)) %>%
                     mutate(respondingTo = ifelse(grepl('YELLEN',respondingTo),'Yellen',
                                                  ifelse(grepl('FED',respondingTo),'Male Fed Chairs',NA))) %>%
                     mutate(respondingTo = relevel(factor(respondingTo),ref = 'Male Fed Chairs')) %>%
                     pull(SENT_combAttack)) - quantile(utterance_level %>%
                                                         filter(all > 30,ind > mind) %>%
                                                         filter(!grepl('FED',opensecretsID)) %>%
                                                         filter(grepl('FED',respondingTo)) %>%
                                                         mutate(respondingTo = ifelse(grepl('YELLEN',respondingTo),'Yellen',
                                                                                      ifelse(grepl('FED',respondingTo),'Male Fed Chairs',NA))) %>%
                                                         mutate(respondingTo = relevel(factor(respondingTo),ref = 'Male Fed Chairs')) %>%
                                                         pull(SENT_combAttack),p = seq(.668897,.668899,by = .00000001))))


# Roughly 67%
toplotRQ <- NULL
for(tau in c(.6688979,seq(.3,.9,by = .1))) {
  # stop()
  mRQ <- rq(SENT_combAttack ~ respondingTo*factor(interruptor) + log(nchars) + log(nchars_lag) + factor(interrupted) + 
              opensecretsID,
            tau = tau,
            data = utterance_level %>%
              filter(all > 30,ind > mind) %>%
              filter(!grepl('FED',opensecretsID)) %>%
              filter(grepl('FED',respondingTo)) %>%
              mutate(respondingTo = ifelse(grepl('YELLEN',respondingTo),'Yellen',
                                           ifelse(grepl('FED',respondingTo),'Male Fed Chairs',NA))) %>%
              mutate(respondingTo = relevel(factor(respondingTo),ref = 'Male Fed Chairs'))
  )
  summary(mRQ)
  tmp <- plot_cme(mRQ,variables = 'respondingTo',condition = 'interruptor',draw = F)
  toplotRQ <- toplotRQ %>%
    bind_rows(tmp %>%
                mutate(tau = tau))
}


pdf('./output/figures/SI_figure_26.pdf',width = 7,height = 5)
toplotRQ %>%
  as_tibble() %>%
  filter(tau != .6688979) %>%
  ggplot(aes(x = tau,y = estimate,color = factor(interruptor),shape = factor(interruptor))) + 
  geom_point(position = position_dodge(width = .005),size = 2.5) + 
  geom_errorbar(aes(ymin = conf.low,ymax = conf.high),
                position = position_dodge(width = .005),width = .005) + 
  geom_hline(yintercept = 0,linetype = 'dashed') + 
  scale_color_manual(name = 'Utterance',values = c('grey50','black'),labels = c('Non-Interruption','Interruption')) + 
  scale_shape_manual(name = 'Utterance',values = c(19,17),labels = c('Non-Interruption','Interruption')) + 
  theme_bw() + 
  theme(legend.position = 'bottom') + 
  labs(x = 'Quantile',
       y = 'Marginal Effect of Yellen\nRelative to Male Fed Chairs',
       title = 'Quantile Regression Results',
       subtitle = 'Marginal effects by decile of aggression outcome')
dev.off()



# SI Figure 27
toAnalBern <- utterance_level %>% 
  filter(nchars > 60,date > as.Date('2006-01-01'),date < as.Date('2018-01-01'))
toAnalPow <- utterance_level %>% 
  filter(nchars > 60,date > as.Date('2014-01-01')) %>%
  mutate(yellenTime = !yellenTime)
toAnalFull <- utterance_level %>% 
  filter(nchars > 60)

dims <- colnames(utterance_level %>% select(matches('SENT_'),-matches('_lag|error')) %>% select(-matches('SEVERE|AUTHOR|LIKELY')))

labels <- NULL
mList <- list()
for(y in c('SENT_combAttack','SENT_combIncoh','SENT_combToxic')) {
  for(p in c('Bern','Pow','Full')) {
    summary(mList[[paste0(y,':',p)]] <- felm(as.formula(paste0('scale(',y,') ~ fedResp*yellenTime + ',
                                                               paste(paste0('topic_',1:100,'_lag'),collapse = ' + '),
                                                               ' + ',
                                                               paste(paste0('scale(',dims[-which(grepl('comb',dims)|dims == y)],'_lag)'),collapse = ' + '),
                                                               ' + poly(scale(log(nchars_lag+1)),3) + poly(scale(log(tot_utterances)),3) + interrupted + poly(year,3)',
                                                               '| opensecretsID + chamber | 0 | opensecretsID')),
                                             get(paste0('toAnal',p))))
    labels <- bind_rows(labels,
                        data.frame(t(summary(mList[[paste0(y,':',p)]])$coefficients['fedResp:yellenTimeTRUE',])) %>%
                          rename(est = Estimate,se = Cluster.s.e.,tstat = t.value,pval = Pr...t..) %>%
                          mutate(dimension = y,
                                 period = p))
  }
}

labels <- labels %>%
  mutate(star = ifelse(pval < .001,'***',
                       ifelse(pval < .01,'**',
                              ifelse(pval < .05,'*','')))) %>%
  mutate(lab = paste0(round(est,3),star,'\n(',round(se,3),')'))

pDescAgg <- toAnalFull %>%
  filter(nchars > 60,
         date > as.Date('2006-01-01')) %>%
  mutate(period = ifelse(date < as.Date('2014-01-01'),'preYellen',
                         ifelse(date < as.Date('2018-01-01'),'PostYellen','Yellen'))) %>%
  group_by(date,fedResp,period) %>%
  summarise(combAttack = mean(SENT_combAttack),
            combIncoh = mean(SENT_combIncoh),
            combToxic = mean(SENT_combToxic),.groups = 'drop') %>%
  mutate(grp = paste0(period,fedResp)) %>%
  ggplot(aes(x = date,y = combAttack,color = factor(fedResp),group = grp)) + 
  geom_vline(xintercept = as.Date('2014-01-01')) +
  geom_vline(xintercept = as.Date('2018-01-01')) + 
  geom_point() + 
  geom_smooth(method = 'lm') + 
  scale_color_manual(name = 'Responding To',values = c('grey50','grey10'),labels = c('Legislators','Fed Chair')) + 
  theme_ridges() + xlab('Date') + ylab('Pr(Attacking)') + 
  theme(legend.position = 'bottom') + 
  ylim(c(0,1)) + 
  ggtitle('DiD: Aggression Index') + 
  annotate(geom = 'label',x = c(as.Date('2014-01-01'),
                                as.Date('2018-01-01')),
           y = c(Inf,Inf),
           vjust = 1,
           size = 3,
           fill = 'white',
           label = labels %>% 
             filter(dimension == 'SENT_combAttack',period %in% c('Bern','Pow')) %>% .$lab)

pDescInc <- toAnalFull %>%
  filter(nchars > 60,
         date > as.Date('2006-01-01')) %>%
  mutate(period = ifelse(date < as.Date('2014-01-01'),'preYellen',
                         ifelse(date < as.Date('2018-01-01'),'PostYellen','Yellen'))) %>%
  group_by(date,fedResp,period) %>%
  summarise(combAttack = mean(SENT_combAttack),
            combIncoh = mean(SENT_combIncoh),
            combToxic = mean(SENT_combToxic),.groups = 'drop') %>%
  mutate(grp = paste0(period,fedResp)) %>%
  ggplot(aes(x = date,y = combIncoh,color = factor(fedResp),group = grp)) + 
  geom_vline(xintercept = as.Date('2014-01-01')) +
  geom_vline(xintercept = as.Date('2018-01-01')) + 
  geom_point() + 
  geom_smooth(method = 'lm') + 
  scale_color_manual(name = 'Responding To',values = c('grey50','grey10'),labels = c('Legislators','Fed Chair')) + 
  theme_ridges() + xlab('Date') + ylab('Pr(Incoherent)') + 
  theme(legend.position = 'bottom') + 
  ylim(c(0,1)) + 
  ggtitle('DiD: Incoherence Index') + 
  annotate(geom = 'label',x = c(as.Date('2014-01-01'),
                                as.Date('2018-01-01')),
           y = c(Inf,Inf),
           vjust = 1,
           size = 3,
           fill = 'white',
           label = labels %>% 
             filter(dimension == 'SENT_combIncoh',period %in% c('Bern','Pow')) %>% .$lab)

pDescTox <- toAnalFull %>%
  filter(nchars > 60,
         date > as.Date('2006-01-01')) %>%
  mutate(period = ifelse(date < as.Date('2014-01-01'),'preYellen',
                         ifelse(date < as.Date('2018-01-01'),'PostYellen','Yellen'))) %>%
  group_by(date,fedResp,period) %>%
  summarise(combAttack = mean(SENT_combAttack),
            combIncoh = mean(SENT_combIncoh),
            combToxic = mean(SENT_combToxic),.groups = 'drop') %>%
  mutate(grp = paste0(period,fedResp)) %>%
  ggplot(aes(x = date,y = combToxic,color = factor(fedResp),group = grp)) + 
  geom_vline(xintercept = as.Date('2014-01-01')) +
  geom_vline(xintercept = as.Date('2018-01-01')) + 
  geom_point() + 
  geom_smooth(method = 'lm') + 
  scale_color_manual(name = 'Responding To',values = c('grey50','grey10'),labels = c('Legislators','Fed Chair')) + 
  theme_ridges() + xlab('Date') + ylab('Pr(Toxic)') + 
  theme(legend.position = 'bottom') + 
  ylim(c(0,1)) + 
  ggtitle('DiD: Toxicity Index') + 
  annotate(geom = 'label',x = c(as.Date('2014-01-01'),
                                as.Date('2018-01-01')),
           y = c(Inf,Inf),
           vjust = 1,
           size = 3,
           # hjust = c(1,0),
           fill = 'white',
           label = labels %>% 
             filter(dimension == 'SENT_combToxic',period %in% c('Bern','Pow')) %>% .$lab)

# COMB
# Is Yellen more toxic?
dyadToAnalYellen <- utterance_level %>%
  arrange(fullInd) %>%
  mutate_at(vars(matches('lag')),function(x) ifelse(is.na(x),0,x)) %>%
  mutate(votepct_rel = ifelse(is.infinite(votepct_rel),1,votepct_rel),
         yellen = ifelse(grepl('YELLEN',opensecretsID),1,0),
         respondingTo = relevel(factor(respondingTo),ref = 'FEDBERNANKE'),
         opensecretsID = relevel(factor(opensecretsID),ref = 'FEDYELLEN'),
         respGroup = ifelse(grepl('YELLEN',respondingTo),'Yellen',
                            ifelse(grepl('BERNANKE',respondingTo),'Bernanke',
                                   ifelse(grepl('GREENSPAN',respondingTo),'Greenspan',
                                          ifelse(grepl('POWELL',respondingTo),'Powell',lag(party)))))) %>% 
  filter(ind > mind,all > 30)



# Are legislators more toxic toward Yellen?
dyadToAnalLegs <- utterance_level %>%
  arrange(fullInd) %>%
  mutate_at(vars(matches('lag')),function(x) ifelse(is.na(x),0,x)) %>%
  mutate(votepct_rel = ifelse(is.infinite(votepct_rel),1,votepct_rel),
         yellen = ifelse(grepl('YELLEN',respondingTo),1,0),
         respondingTo = relevel(factor(respondingTo),ref = 'FEDYELLEN'),
         opensecretsID = relevel(factor(opensecretsID),ref = 'FEDYELLEN'),
         respGroup = relevel(factor(ifelse(grepl('YELLEN',respondingTo),'Yellen',
                                           ifelse(grepl('BERNANKE',respondingTo),'Bernanke',
                                                  ifelse(grepl('GREENSPAN',respondingTo),'Greenspan',
                                                         ifelse(grepl('POWELL',respondingTo),'Powell',lag(party)))))),ref = 'Yellen')) %>% 
  filter(ind > mind,all > 30)

# Is yellen more toxic, attacking, or incoherent compared to other fed chairs?
toplot <- NULL
mList <- mListRev <- list()
for(y in dims) {
  (mList[[y]] <- feols(as.formula(paste0('scale(',y,') ~ yellen + ',
                                         # paste0(y,'_lag'),
                                         # ' + ',
                                         paste(paste0('topic_',1:100,'_lag'),collapse = ' + '),
                                         ' + ',
                                         paste(paste0('scale(',dims[-which(grepl('comb',dims))],'_lag)'),collapse = ' + '),
                                         ' + poly(log(nchars_lag+1),3) + scale(log(tot_utterances)) + interrupted',
                                         '| respondingTo + chamber')),
                       dyadToAnalYellen %>%
                         filter(grepl('FED',opensecretsID),
                                !grepl('FED',respondingTo)),cluster = c('respondingTo','docID')))
  
  toplot <- bind_rows(toplot,mList[[y]]$coeftable %>%
                        data.frame() %>%
                        rename(est = Estimate,se = Std..Error,tstat = t.value,pval = Pr...t..) %>%
                        slice(1) %>%
                        mutate(out = y,
                               model = 'fedTone') %>% as_tibble())
  
  (mListRev[[y]] <- feols(as.formula(paste0('scale(',y,') ~ yellen + ',
                                            # paste0(y,'_lag'),
                                            # ' + ',
                                            paste(paste0('topic_',1:100,'_lag'),collapse = ' + '),
                                            ' + ',
                                            paste(paste0('scale(',dims[-which(grepl('comb',dims))],'_lag)'),collapse = ' + '),
                                            ' + poly(log(nchars_lag+1),3) + scale(log(tot_utterances)) + interrupted',
                                            '| opensecretsID + chamber')),
                          dyadToAnalLegs %>%
                            filter(grepl('FED',respondingTo),
                                   !grepl('FED',opensecretsID)),cluster = c('respondingTo','docID')))
  
  toplot <- bind_rows(toplot,mListRev[[y]]$coeftable %>%
                        data.frame() %>%
                        rename(est = Estimate,se = Std..Error,tstat = t.value,pval = Pr...t..) %>%
                        slice(1) %>%
                        mutate(out = y,
                               model = 'legTone') %>% as_tibble())
}

pComb <- toplot %>%
  filter(model == 'legTone') %>%
  mutate(sentType = ifelse(out %in% paste0('SENT_',c('ATTACK_ON_AUTHOR','ATTACK_ON_COMMENTER','IDENTITY_ATTACK','INSULT','THREAT','combAttack')),'Aggression',
                           ifelse(out %in% paste0('SENT_',c('TOXICITY','SEVERE_TOXICITY','PROFANITY','SEXUALLY_EXPLICIT','FLIRTATION','OBSCENE','INFLAMMATORY','combToxic')),'Toxicity',
                                  ifelse(out %in% paste0('SENT_',c('INCOHERENT','UNSUBSTANTIAL','combIncoh')),'Incoherence',NA))),
         out = gsub('Attack On Commenter','Attack',
                    gsub('Combattack','INDEX: Aggression',
                         gsub('Combincoh','INDEX: Incoherence',
                              gsub('Combtoxic','INDEX: Toxicity',str_to_title(gsub('_',' ',gsub('SENT_','',out)))))))) %>%
  mutate(index = ifelse(grepl('INDEX',out),'index','dimension')) %>%
  filter(!is.na(sentType),
         !out %in% c('Attack On Author','Severe Toxicity'),
         !grepl('INDEX:',out)) %>%
  ggplot(aes(x = est,y = reorder(out,est),shape = index,fill = index,size = index)) + 
  geom_errorbarh(aes(xmin = est - 2*se,xmax = est + 2*se),height = .1,size = .5) + 
  geom_point() + 
  geom_vline(xintercept = 0,linetype = 'dashed') + 
  theme_ridges() + 
  scale_shape_manual(values = c(19,21)) + 
  scale_size_manual(values = c(2,4)) + 
  scale_fill_manual(values = c('white','white')) + 
  facet_grid(sentType~.,scales = 'free',shrink = TRUE,space = 'free') + 
  theme(legend.position = 'none') + 
  xlab("Yellen coefficient") + 
  ylab('') + 
  labs(title = 'Interlocutor FE',subtitle = 'Yellen Relative to Male Chairs')

leg <- get_legend(pDescAgg + theme(legend.text = element_text(size = 8),legend.title = element_text(size = 9)))
pdf('./output/figures/SI_figure_27.pdf',width = 7,height = 7)
ggdraw() + 
  draw_plot(pDescAgg + theme(legend.position = 'none',
                             plot.margin = unit(c(0,0,-.5,0), "cm")) + xlab(''),x = 0,y = .68,width = .4,height = .32) + 
  draw_plot(pDescInc + theme(legend.position = 'none',
                             plot.margin = unit(c(0,0,-.5,0), "cm")) + xlab(''),x = 0,y = .36,width = .4,height = .32) + 
  draw_plot(pDescTox + theme(legend.position = 'none',
                             plot.margin = unit(c(0,0,-.5,0), "cm")) + xlab(''),x = 0,y = .04,width = .4,height = .32) + 
  draw_plot(leg,x = .05,y = 0,width = .4,height = .05) + 
  draw_plot(pComb + theme(legend.position = 'none',
                          plot.margin = unit(c(0,0,0,0), "cm")),x = .4,y = 0,width = .6,height = 1)
dev.off()


# SI Figure 28
# Daughters
toplot <- NULL
mList <- list()
for(dim in colnames(utterance_level %>% select(matches('SENT_',ignore.case = F)))) {
  cat(dim,'\n')
  if(grepl('error|_lag|DEM|GOP',dim)) { next }
  mList[[dim]] <- feols(as.formula(paste0(dim,' ~ respondingTo*anyDaughters + ',
                                          paste(paste0('topic_',1:100,'_lag'),collapse = ' + '),
                                          ' + ',
                                          paste(paste0('scale(',dims[-which(grepl('comb',dims)|dims == y)],'_lag)'),collapse = ' + '),
                                          ' + poly(log(nchars_lag+1),3) + poly(log(tot_utterances),3) + interrupted + poly(year,3)',
                                          '| opensecretsID + chamber')),
                        utterance_level %>%
                          filter(grepl("FED",respondingTo),
                                 nchars > 60) %>%
                          mutate(yellen = ifelse(grepl('YELLEN',respondingTo),1,0),
                                 anyDaughters = ifelse(nDaughters > 0,1,0),
                                 propDaughters = nDaughters / nKids,
                                 Male = ifelse(gender == 'M',1,0)))
  
  tmp <- plot_cme(mList[[dim]],variables = 'respondingTo',by = 'anyDaughters',draw = F)
  
  toplot <- toplot %>% 
    bind_rows(tmp %>%
                mutate(dim = dim))
  
}

toadd <- NULL
for(y in names(mList)) {
  toadd <- toadd %>%
    bind_rows(summary(mList[[y]])$coeftable %>%
                data.frame() %>%
                mutate(covs = row.names(.)) %>%
                as_tibble() %>%
                filter(grepl(':',covs)) %>%
                rename(est = Estimate,se = Std..Error,tstat = t.value,pval = Pr...t..) %>%
                mutate(contrast = str_to_title(gsub('respondingToFED|:anyDaughters','',covs))) %>%
                mutate(dim = y))
}

toplot2 <- toplot %>%
  as_tibble() %>%
  filter(grepl('comb',dim)) %>%
  mutate(contrast = str_to_title(gsub('mean\\(|\\)| - |FEDBERNANKE|FED','',contrast))) %>%
  left_join(toadd) %>%
  mutate(dim = ifelse(grepl('Attack',dim),'Aggression',
                      ifelse(grepl('Incoh',dim),'Incoherent','Toxic'))) %>%
  mutate(stars = ifelse(pval < .001,'***',
                        ifelse(pval < .01,'**',
                               ifelse(pval < .05,'*',
                                      ifelse(pval < .1,'.','')))))

pdf('./output/figures/SI_figure_28.pdf',width = 7,height = 5)
toplot2 %>%
  ggplot(aes(x = factor(anyDaughters),
             y = estimate,group = dim)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = conf.low,ymax = conf.high),width = 0) + 
  geom_line() + facet_grid(dim~contrast) + 
  scale_x_discrete(breaks = c(0,1),
                   labels = c('None','Any')) + 
  geom_hline(yintercept = 0,linetype = 'dashed') + 
  geom_label(data = toplot2 %>%
               filter(anyDaughters == 1),
             x = 1.5,y = .05,aes(label = paste0(round(est,3),stars,'\n(',round(se,3),')'))) + 
  theme_bw() + 
  labs(x = 'Daughters',
       y = 'Effect of Chair\nCompared to Bernanke')
dev.off()

# EOF